*This program runs a probit and calculates the associated marginal effects
*for stocks, own business and home in the 11 SHARE countries
*It's used on the paper on wealth counterfactuals 
*for the final verification for the Review of Economics and Statistics
*June 2011 - Dimitris Christelis

version 9.2
cap clear all
macro drop _all
clear
set mem 500m
set more off

*Paths
do "c:/foldp/foldp.do"
global projf "${foldp}/Wealth Counterfactuals/Verification"

*Tolerance
sca tolr=1e-6

*Weights
global wgt wgth

*ATTENTION: The base country (Germany in our case) has to be put first in the sequence
global eureg1 de se dk nl be fr ch at it es gr
global eureg2 Germany Sweden Denmark Netherlands Belgium France Switzerland Austria ///
   Italy Spain Greece
*Country numbers: they have to follow exactly the same sequence as the country names
global eureg3 3 1 2 4 5 6 7 8 9 10 11   
local nas: word count $eureg1

global assetl1 rfin home ownb secdebt liab tdebt
global assetl2 hrfin hhome hownb hsecdebt hliab htdebt
local nar: word count $assetl1

*The loop for the asset must come before the loop for the countries
 
forvalues ll=1/`nar' {

   *Globals of asset
   global as: word `ll' of $assetl1
   global hhas: word `ll' of $assetl2

   forvalues kk=1/`nas' {

   *Globals of country
   global re: word `kk' of $eureg1
   global fnre: word `kk' of $eureg2
   global cn: word `kk' of $eureg3

   *Logging
   cap log close
   log using "${projf}/Programs/SHARE/${re}_${as}.log", replace

   global pdate "$S_DATE" 
   global ptime "$S_TIME"
   display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"

   *Opening HRS data to generate quartiles
   use "${projf}/Data/HRS/HRS_DG.dta", clear

   *Checking the sample weights
   assert abs(wgth-wgtach)<tolr
   drop wgtach

   *Keeping only heads
   keep if head==1 & wgth>0 & wgth!=.

   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov /*+ hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Creating wealth quartiles 
   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_us = r(p25)
   sca p50_w_us = r(p50)
   sca p75_w_us = r(p75)

   *Creating income quartiles
   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_us = r(p25)
   sca p50_i_us = r(p50)
   sca p75_i_us = r(p75)


   *Reading the parameters, variable lists etc.
   do "${projf}/Programs/Parameters.do"

   *Reading the SHARE data
   use "${projf}/Data/SHARE/SHARE_DG.dta", clear

   *Keeping only the country of interest
   keep if country==${cn} & ${wgt}>0 & ${wgt}!=.

   *Keeping only heads
   keep if head==1

   *Total debt
   qui gen double htdebtv = hmortv+hliabv /* +hhmlov*/
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt
   qui ge double hsecdebtv = hmortv
   qui ge double hsecdebto = (hsecdebtv>tolr)

  
   *Rescaling age
   qui replace age=age/sc_age
   qui gen double agesq=age^2

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Bad health dummy
   qui gen byte badhs=(srhealtha>=4)

   *College dummy
   qui gen byte college=(edu==3)
   
   *High school dummy
   qui gen byte highs=(edu==2)

   *Couple dummy
   qui gen byte couple=(mstatus==1)
   
   *Widow dummy
   qui gen byte widow=(mstatus==3)
   
   *Never married
   qui gen byte nevmar=(mstatus==4)
      
   *Creating wealth quartiles for financial and real assets, both with respect
   *to US and DE
   if ${cn}==3 {
      qui su nw if head==1 [aw=wgth], detail
      sca p25_w_de = r(p25)
      sca p50_w_de = r(p50)
      sca p75_w_de = r(p75)
   }
   qui gen double wqnw1_qus=(nw<=scalar(p25_w_us))
   qui gen double wqnw2_qus=(nw>scalar(p25_w_us) & nw<=scalar(p50_w_us))
   qui gen double wqnw3_qus=(nw>scalar(p50_w_us) & nw<=scalar(p75_w_us))
   qui gen double wqnw4_qus=(nw>scalar(p75_w_us))
   qui gen double wqnw1_qde=(nw<=scalar(p25_w_de))
   qui gen double wqnw2_qde=(nw>scalar(p25_w_de) & nw<=scalar(p50_w_de))
   qui gen double wqnw3_qde=(nw>scalar(p50_w_de) & nw<=scalar(p75_w_de))
   qui gen double wqnw4_qde=(nw>scalar(p75_w_de))
   
   if ${cn}==3 {   
      qui su inc if head==1 [aw=wgth], detail
      sca p25_i_de = r(p25)
      sca p50_i_de = r(p50)
      sca p75_i_de = r(p75)
   }
   qui gen double wqinc1_qus=(inc<=scalar(p25_i_us))
   qui gen double wqinc2_qus=(inc>scalar(p25_i_us) & inc<=scalar(p50_i_us))
   qui gen double wqinc3_qus=(inc>scalar(p50_i_us) & inc<=scalar(p75_i_us))
   qui gen double wqinc4_qus=(inc>scalar(p75_i_us))
   qui gen double wqinc1_qde=(inc<=scalar(p25_i_de))
   qui gen double wqinc2_qde=(inc>scalar(p25_i_de) & inc<=scalar(p50_i_de))
   qui gen double wqinc3_qde=(inc>scalar(p50_i_de) & inc<=scalar(p75_i_de))
   qui gen double wqinc4_qde=(inc>scalar(p75_i_de))


   *Saving the dataset
   qui save "${projf}/Results/SHARE/${re}/${as}/temp.dta", replace

   *Running 2 regressions, using both definitions of quartiles
   foreach cb in qus qde {
   
      use "${projf}/Results/SHARE/${re}/${as}/temp.dta", clear
      probit ${hhas}o ${basiccovs_`cb'}, r

      *Coefficient matrix
      mat olco=e(b)
  
      *Augmenting the coefficient matrix by an empty column, the number of obs and the log likelihood
      mat blank = .
      mat colnames blank=blank
  
      mat olN = e(N)
      mat colnames olN=nobs
  
      mat ill = e(ll)
      mat colnames ill=loglik

      mat olco=[olco,blank,olN,ill]
       
      *Var-Covar matrix
      mat olvar=e(V)

      *Saving matrix of regression results as a text file
  
      *Number of variables
      local nc=colsof(olvar)
      global nc=`nc'
      matrix olout=J(`nc'+3,6,.)

      *Checking that the number of regressors actually used is the same 
      *for all countries for a given asset (i.e. no variable is dropped)
      if "${re}"=="de" {
         global nc_de_${as}_`cb'=`nc'
      }
      if "${re}"~="de" {  
         assert `nc'==${nc_de_${as}_`cb'}
      }

  
      *Putting in the regression coefficients 
      forvalues k=1/`nc' {
         matrix olout[`k',1]=olco[1,`k']
         matrix olout[`k',2]=sqrt(olvar[`k',`k'])
         matrix olout[`k',3]=olco[1,`k']/sqrt(olvar[`k',`k'])
         matrix olout[`k',4]=2*ttail(e(N) - `nc',abs(olout[`k',3]))
         matrix olout[`k',5] = olout[`k',1] - invttail(e(N) - `nc',sigl/2)*olout[`k',2]
         matrix olout[`k',6] = olout[`k',1] + invttail(e(N) - `nc',sigl/2)*olout[`k',2]
      }
  
      *Adding number of observations used in first stage probit
      matrix olout [`nc'+2,1]=e(N)
      *Adding log likelihood
      matrix olout [`nc'+3,1]=e(ll)

      *Labels of matrix of regression output
      *Column labels
      local colmat coeff se tstat pvalue ci_low ci_high
      matrix colnames olout = `colmat' 
      *Row labels
      matrix rownames olout = ${basiccovs_`cb'} _cons blank nobs loglik
      *Displaying and saving the matrix of results
      mat list olout,format(%15.5f)
      mat2txt, mat(olout) saving("${projf}/Results/SHARE/${re}/${as}/probit_est_${re}_${as}_`cb'.txt") replace format(%15.5f) ///
        title("Regression for ${as} in ${fnre}, using `cb' quartiles")
  
      preserve
      *Saving matrices of coefficients and correlations and their covariance matrices
      *for each implicate
      drop _all
      svmat double olco, names(eqcol)
      qui save "${projf}/Results/SHARE/${re}/${as}/probit_coeff_${re}_${as}_`cb'.dta", replace
      drop _all
      svmat double olvar, names(eqcol)
      qui save "${projf}/Results/SHARE/${re}/${as}/probit_var_${re}_${as}_`cb'.dta", replace
      restore  

      *Bootstrap program

      *Creating a matrix with number of columsn equal to 1 + # bootstrap runs, storing
      *the main sample estimates in the first line, together with the proportion of
      *hhds owning the asset
      qui su ${hhas}o if head==1 [aw=$wgt], meanonly
      sca ${as}p=r(mean)
      mat olco=olco[1,1..${nc}]
      mat ${as}coeffall=(olco,${as}p)
      local nch=$nc+1
      local bn ""
      forvalues lll=1/`nch' {
        local bn `bn' b`lll'
      }
      mat colnames ${as}coeffall = `bn'   
      qui save "${projf}/Results/SHARE/${re}/${as}/temp2.dta", replace
  
      *Subsequent runs, beginning of the bootstrap loop
      forvalues i=1/$finsim {

         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
           
         use "${projf}/Results/SHARE/${re}/${as}/temp2.dta", clear
         *Fixing the seed before each resampling
         local sd=1234+30*`i'+`kk'
         set seed `sd'
         bsample
         
         *Probit regression in the bootstrapped sample
         if `i'==1 {
            probit ${hhas}o ${basiccovs_`cb'}, r
         }
         if `i'>1 {
            qui probit ${hhas}o ${basiccovs_`cb'}, r
         }

         mat a=e(V)
         local nc`i'=colsof(a)
               
         *If no variables are dropped from the probit then we proceed with updating the matrix
         if `nc`i''==$nc {
            qui su ${hhas}o if head==1 [aw=$wgt], meanonly
            sca ${as}p=r(mean)
            mat olco=e(b)
            mat ${as}coeff`sr'=(olco,${as}p)
            mat ${as}coeffall=(${as}coeffall \ ${as}coeff`sr')

            forvalues lll=1/10 {
               if `sr'==$nsim*(`lll'/10) {
                  noisily di in yellow `sr' in yellow "*" _continue 
               }
            }

            mat drop ${as}coeff`sr'
         }
         
         *If some variables are dropped from the probit in the bootstrapped samples 
         *we put back the index to its previous value
         if `nc`i''~=$nc {
            local sr=`sr'-1
         }
               
         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         if `sr'==$nsim {
            di in yellow "End of ${nsim} valid boostrap runs, total runs equal to `i'"
            continue, break
         }
         
      }  /* End of loop over bootstrap runs */

      *Saving the matrix of bootstrap coefficients and sample proportions
      drop _all
      svmat double ${as}coeffall, names(eqcol)
      qui save "${projf}/Results/SHARE/${re}/${as}/boot_drawcoeff_${re}_${as}_`cb'.dta", replace

        
      ***********************************************************
      *Calculation of marginal effects and associated magnitudes*
      ***********************************************************

      *Reading the matrices saved from the estimation part, to retrieve the # of obs
      use "${projf}/Results/SHARE/${re}/${as}/probit_coeff_${re}_${as}_`cb'.dta", clear
      mkmat _all , mat(a) 
      local nc=colsof(a)
      global nc=`nc'-3
      sca nobs_`cb'=a[1,${nc}+2]


      *Reading the dataset
      use "${projf}/Results/SHARE/${re}/${as}/temp.dta", clear

      *Incrementing the continuous variables by the appropriate magnitudes
      *to compute marginal effects and elasticities
      *Variables that are not used in elasticity form
      foreach y in $conlist_meff {
         qui gen double `y'_pl=`y'+(ch_`y'/sc_`y')
      }   

      merge using "${projf}/Results/SHARE/${re}/${as}/boot_drawcoeff_${re}_${as}_`cb'.dta"   
      drop _merge   

      *Putting the values of the drawn coefficients equal to the dummy value from missing
      *in rows lower than the number of bootstraps, so that there are no missing values
      local nch = $nc + 1
      forvalues h=1/`nch' {   
         qui replace _b`h'=$dumval if _n>${nsim}+1  
      }      
   
      *Renaming the coefficients  
      local ncl = $nc - 1
      forvalues h=1/`ncl' {   
         local g: word `h' of ${basiccovs_`cb'}
         rename _b`h' b_`g'_`cb'   
      }   

      *Renaming the coefficient of the constant 
      *set trace on
      rename _b${nc} b_cons_`cb'
   
      local nch = $nc + 1
      rename _b`nch' pr_${re}_${as}_`cb'_own
   
      *Merging with the draws from the US
      if "`cb'"=="qus" {
         merge using "${projf}/Results/HRS/Country/${as}/boot_drawcoeff_usall_${as}.dta"   
      	 drop _merge   

      	 *Making sure that the number of coefficients in the SHARE regression is the 
      	 *same as the one in the US regression
      	 preserve
      	 use "${projf}/Results/HRS/Country/${as}/probit_var_usall_${as}.dta", clear
      	 mkmat _all , mat(olvar) 
      	 local chnc=colsof(olvar)
      	 assert `chnc'==${nc_de_${as}_`cb'}
      	 restore

      	 *Global of number of variables
      	 global nc=colsof(olvar)


      	 *Putting the values of the drawn coefficients equal to the dummy value from missing
      	 *in rows lower than the number of bootstraps, so that there are no missing values
      	 local nch = $nc + 1
      	 forvalues h=1/`nch' {   
      	    qui replace _b`h'=$dumval if _n>${nsim}+1  
      	 }      
   
      	 *Renaming the coefficients  
      	 local ncl = $nc - 1
      	 forvalues h=1/`ncl' {   
      	    local g: word `h' of ${basiccovs_qus}
      	    rename _b`h' b_usall_`g'   
      	 }   

      	 *Renaming the coefficient of the constant 
      	 *set trace on
      	 rename _b${nc} b_usall_cons
   
      	 local nch = $nc + 1
      	 rename _b`nch' pr_usall_${as}_own
      }
      
      *Merging with the coefficients and the probabilities from Germany
      if "`cb'"=="qde" {
         if "${re}"~="de" {

      	    *Coefficients
      	    merge using "${projf}/Results/SHARE/de/${as}/boot_drawcoeff_de_${as}_qde.dta"   
      	    drop _merge   

      	    *Putting the values of the drawn coefficients equal to the dummy value from missing
      	    *in rows lower than the number of bootstraps, so that there are no missing values
      	    local nch = $nc + 1
      	    forvalues h=1/`nch' {   
      	       qui replace _b`h'=$dumval if _n>${nsim}+1  
      	    }      
   
      	    *Renaming the coefficients  
      	    local ncl = $nc - 1
      	    forvalues h=1/`ncl' {   
      	       local g: word `h' of ${basiccovs_qde}
      	       rename _b`h' b_de_`g'_qde   
      	    }   

      	    *Renaming the coefficient of the constant 
      	    *set trace on
      	    rename _b${nc} b_de_cons_qde
   
      	    local nch = $nc + 1
      	    rename _b`nch' pr_de_${as}_qde_own
      	 }         
      }
      
      *Loop over bootstrap coefficients
      forvalues i=1/$finsim { 

         *Defining an index of valid bootstrap runs that will determine whether the 
         *loop finishes or not
         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
         *If there were any condition on the coefficients that was violated by the current
         *draw then the index of valid bootstrap runs would be decreased by one at this point
         *in the program. For the pooled probit there is no such condition and thus the
         *index of valid bootstrap runs (sr) is equal to the on of total bootstrap runs (i)
   
         *Creating the linear prediction for the current draw of the coefficients
         qui gen double xb0_`cb' = b_cons_`cb'[`i']
         local lll ${basiccovs_`cb'}
         foreach y in `lll' {
            qui replace xb0_`cb' = xb0_`cb' + b_`y'_`cb'[`i']*`y' 
         }
      
         *Creating a counterfactual prediction using the coefficients from the US
         if "`cb'"=="qus" {
            qui gen double xb0_usall = b_usall_cons[`i']
            foreach y in $basiccovs_qus {
               qui replace xb0_usall = xb0_usall + b_usall_`y'[`i']*`y' 
            }
	 }

         *Creating a counterfactual prediction using the coefficients from Germany
         if "`cb'"=="qde" {
            if "${re}"~="de" {         
               qui gen double xb0_de_qde = b_de_cons_qde[`i']
               foreach y in $basiccovs_qde {
                  qui replace xb0_de_qde = xb0_de_qde + b_de_`y'_qde[`i']*`y' 
               }
            }            
         }
      
         *******************************************************************
         *Calculating the average of probabilities across individual obs
         if "`cb'"=="qus" {   
            qui gen double dp_usall=normal(xb0_usall)
            qui su dp_usall [aw=$wgt], meanonly
            sca r0m=r(mean)
         }
         
         *Probability using Germany coefficients
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               qui gen double dp_de_qde=normal(xb0_de_qde)
               qui su dp_de_qde [aw=$wgt], meanonly
               sca r2m_qde=r(mean)
            }
         }

         if `i'==1 {
             
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients
               qui gen double pr_${re}_${as}_qus_cus=r0m
             
               *Difference between US ownership proportion and
               *own ownership proportion
               qui gen double dpus_${re}_${as}_qus_tot=pr_usall_${as}_own - ///
                  pr_${re}_${as}_qus_own
             
               *Difference between US and own probability calculated using
               *US coefficients (characteristics effect)
               qui gen double dpus_${re}_${as}_qus_char=pr_usall_${as}_own[1] - ///   
                  pr_${re}_${as}_qus_cus[1]
             
               *Difference between own probability calculated using
               *US coefficients and own probability (coefficient effect)
               qui gen double dpus_${re}_${as}_qus_coeff=pr_${re}_${as}_qus_cus[1] - ///   
                  pr_${re}_${as}_qus_own[1]
            } 

            if "`cb'"=="qde" {
               if "${re}"~="de" {
            
             	  *Counterfactual probability using German coefficients
             	  qui gen double pr_${re}_${as}_qde_cde=r2m_qde
             
             	  *Difference between mean German predicted probability and
             	  *mean own predicted probability
             	  qui gen double dpde_${re}_${as}_qde_tot=pr_de_${as}_qde_own - ///
             	     pr_${re}_${as}_qde_own
             
             	  *Difference between German and own probability calculated using
             	  *German coefficients (characteristics effect)
             	  qui gen double dpde_${re}_${as}_qde_char = pr_de_${as}_qde_own[1] - ///
             	     pr_${re}_${as}_qde_cde[1]  
             
             	  *Difference between own probability calculated using
             	  *German coefficients and own probability (coefficient effect)
             	  qui gen double dpde_${re}_${as}_qde_coeff = pr_${re}_${as}_qde_cde[1] - ///
             	    pr_${re}_${as}_qde_own[1]
               }
            }
         }
         
         *Putting the probability at the line indexed by
         *the valid bootstrap run (sr), and not the general
         *run (i)
         
         if `i'>1 {
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients   
               qui replace pr_${re}_${as}_qus_cus=r0m in `sr'   
               assert pr_${re}_${as}_qus_cus<$dumval if _n==`sr'   

               *Difference between US and own probability calculated using   
               *US coefficients (characteristics effect)   
               qui replace dpus_${re}_${as}_qus_char=pr_usall_${as}_own[`sr'] - ///   
                  pr_${re}_${as}_qus_cus[`sr'] in `sr'   

               *Difference between own probability calculated using   
               *US coefficients and own probability (coefficient effect)   
               qui replace dpus_${re}_${as}_qus_coeff=pr_${re}_${as}_qus_cus[`sr'] - ///   
                  pr_${re}_${as}_qus_own[`sr'] in `sr'   
            }

            if "`cb'"=="qde" {
               if "${re}"~="de" {
                  *Counterfactual probability using German coefficients
               	  qui replace pr_${re}_${as}_qde_cde=r2m_qde in `sr'
               	  assert pr_${re}_${as}_qde_cde<$dumval if _n==`sr'
         
               	  *Difference between German and own probability calculated using
               	  *German coefficients (characteristics effect)
               	  qui replace dpde_${re}_${as}_qde_char=pr_de_${as}_qde_own[`sr'] - ///
               	     pr_${re}_${as}_qde_cde[`sr'] in `sr'

               	  *Difference between own probability calculated using
               	  *German coefficients and own probability (coefficient effect)
               	  qui replace dpde_${re}_${as}_qde_coeff=pr_${re}_${as}_qde_cde[`sr'] - ///
               	     pr_${re}_${as}_qde_own[`sr'] in `sr'
               } 
            }
         }
         
         if "`cb'"=="qus" {
            drop dp_usall xb0_usall
         }
         
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               drop dp_de_qde xb0_de_qde
            }   
         }   
         
         ******************************************************************************************        	 
         ************* 	
         *0/1 Dummies* 	
         ************* 	
    	
         foreach y in $dumlist_meff /* $dumlist */ {           	
            *Putting all dummies to zero          	
            *Removing the effect of the dummy.         	
            qui gen double xb_c1_`cb' = xb0_`cb' - b_`y'_`cb'[`i']*`y'    	        	
            assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
            *Putting back the dummy equal to 1 	
            qui gen double xb_c2_`cb' = xb_c1_`cb' + b_`y'_`cb'[`i']    	        	
            assert xb_c2_`cb'~=. if xb0_`cb'~=.   	        	
    	
            *Loop over the marginal effect and percentage variation 	
            forvalues tt=1/$wcp { 	
               local wpf: word `tt' of ${prfx} 	
     	
      	       if `tt'==1 {    	
                  qui gen double pdif_`y'_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
               }    	
    	       if `tt'==2 {             	
	          qui ge double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	             normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
	       }    	
	       assert pdif_`y'_`cb'~=. if xb0_`cb'~=.   	
               qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                    	
   	       if `i'==1 {                   	
	          qui ge double `wpf'_`y'_`cb' = $dumval                   	
	       }                   	
               *Putting the effect at the appropriate line in the file
               qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                      	
	       assert `wpf'_`y'_`cb'<$dumval if _n==`sr' 	
               drop pdif_`y'_`cb'
	     		         	
            }   /* End of loop over magnitudes */ 	
            drop xb_c1_`cb' xb_c2_`cb' 	
         }  /* End of loop over list of 0/1 dummies */        	
	 
         ******************************************************************************************        	 
         ************************ 	
         *Interdependent Dummies* 	
         ************************ 	
         foreach g in ${groupno_meff_`cb'} /* forvalues g=1/$ng */ {  

            *Putting all dummies in a group equal to zero            
            qui gen double xb_c1_`cb' = xb0_`cb'
            local m ${group`g'}
            foreach y in `m'  {            
               qui replace xb_c1_`cb' = xb_c1_`cb' - b_`y'_`cb'[`i']*`y'            
            }	
            
            local m ${group`g'}
            foreach y in `m'  {     
               *Putting back the dummy equal to 1
     	       qui gen double xb_c2_`cb'= xb_c1_`cb' + b_`y'_`cb'[`i'] 

               *Loop over the elasticity and the 4 intervals
               forvalues tt=1/$wcp {
                  local wpf: word `tt' of ${prfx}    
          
	          if `tt'==1 {   
                     qui gen double pdif_`y'_`cb' = normal(xb_c2_`cb') - normal(xb_c1_`cb')     
	          }   
    	          if `tt'==2 {   
	             qui gen double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	                normal(xb_c1_`cb'))/normal(xb_c1_`cb')                     
	          }   
  	          assert pdif_`y'_`cb'~=. if xb0_`cb'~=.

  	          qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                     
  	          if `i'==1 {                     
	             qui ge double `wpf'_`y'_`cb' = $dumval                     
	          }                     
                  qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                     
	          assert `wpf'_`y'_`cb'<$dumval if _n==`sr'   
                  drop pdif_`y'_`cb' 
  
              }  /* End of loop over magnitudes */
               drop xb_c2_`cb' 
            }  /* End of loop over variables of a given group */       
            drop xb_c1_`cb'
         }  /* End of loop over groups */       

         ******************************************************************************************        	 
         ********************** 	
         *Continuous variables* 	
         ********************** 	
    	
         foreach y in $conlist_meff {           	
            qui gen double xb_c1_`cb' = xb0_`cb' 
            assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
            *Incrementing by the already determined sum	
            qui gen double xb_c2_`cb' = xb_c1_`cb' + b_`y'_`cb'[`i']*(`y'_pl - `y')    	        	
            assert xb_c2~=. if xb0_`cb'~=.   	        	
    	
            *Loop over the marginal effect and percentage variation 	
            forvalues tt=1/$wcp { 	
               local wpf: word `tt' of ${prfx} 	
     	
      	       if `tt'==1 {    	
                    qui gen double pdif_`y'_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
               }    	
    	       if `tt'==2 {             	
	          qui ge double pdif_`y'_`cb' = (normal(xb_c2_`cb') - ///
	             normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
	       }    	
	       assert pdif_`y'_`cb'~=. if xb0_`cb'~=.   	
               qui sum pdif_`y'_`cb' if xb0_`cb'~=. [aw = $wgt], meanonly                    	
   	       if `i'==1 {                   	
	          qui ge double `wpf'_`y'_`cb' = $dumval                   	
	       }                   	
               *Putting the effect at the appropriate line in the file
               qui replace `wpf'_`y'_`cb' = r(mean) in `sr'                      	
	       assert `wpf'_`y'_`cb'<$dumval if _n==`sr' 	
  	       drop pdif_`y'_`cb'	         	

            }   /* End of loop over magnitudes */ 	
            drop xb_c1_`cb' xb_c2_`cb' 	
         }  /* End of loop over list of continuous variables */        	
         
         ******************************************************************************************        	 
         *****
         *Age* 	
         *****
         qui gen double xb_c1_`cb' = xb0_`cb' 
         assert xb_c1_`cb'~=. if xb0_`cb'~=.   	        	
	    	        	
         *Incrementing by the already determined sum	
         qui gen double xb_c2_`cb' = xb_c1_`cb' + b_age_`cb'[`i']*(ch_age/sc_age) + ///
            b_agesq_`cb'[`i']*(2*age*(ch_age/sc_age) + ((ch_age/sc_age)^2))
         assert xb_c2_`cb'~=. if xb0_`cb'~=.   	        	
    	
         *Loop over the marginal effect and percentage variation 	
         forvalues tt=1/$wcp { 	
            local wpf: word `tt' of ${prfx} 	
     	
            if `tt'==1 {    	
                 qui gen double pdif_age_`cb'= normal(xb_c2_`cb') - normal(xb_c1_`cb') 	
            }    	
            if `tt'==2 {             	
               qui ge double pdif_age_`cb' = (normal(xb_c2_`cb') - ///
                 normal(xb_c1_`cb'))/normal(xb_c1_`cb')           	  
            }    	
            assert pdif_age_`cb'~=.    	
            qui sum pdif_age_`cb' [aw = $wgt], meanonly                    	
            if `i'==1 {                   	
               qui ge double `wpf'_age_`cb' = $dumval                   	
            }                   	
            *Putting the effect at the appropriate line in the file
            qui replace `wpf'_age_`cb' = r(mean) in `sr'                      	
	    assert `wpf'_age_`cb'<$dumval if _n==`sr' 	
            drop pdif_age_`cb'	     		         	

         }   /* End of loop over magnitudes */ 	
         drop xb_c1_`cb' xb_c2_`cb' 	

         **********************************************************************************	 
         *Dropping linear prediction based on current draw
         drop xb0_`cb' 
      
         forvalues lll=1/10 {
            if `sr'==$nsim*(`lll'/10) {
               di in yellow "Iteration `sr', ${as}"
            }
         }

         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         local ncpp=$nsim+1
         if `sr'==`ncpp' {
            di in yellow "End of ${nsim} valid bootstrap runs, total runs equal to `i'"
            continue, break
         }
      
      }  /* End of loop over bootstrap runs */

      *Dropping obs beyond the line of the last valid bootstrap run
      drop if _n>$nsim+1
   
      *Global of marginal effects, putting first the probabilities of interest
      if "`cb'"=="qus" {
         local me pr_${re}_${as}_qus_own pr_${re}_${as}_qus_cus dpus_${re}_${as}_qus_tot ///
            dpus_${re}_${as}_qus_coeff dpus_${re}_${as}_qus_char
      }
      if "`cb'"=="qde" {
         if "${re}"=="de" {
            local me pr_${re}_${as}_qde_own
         } 
         if "${re}"~="de" {
            local me pr_${re}_${as}_qde_own pr_${re}_${as}_qde_cde dpde_${re}_${as}_qde_tot ///
               dpde_${re}_${as}_qde_coeff dpde_${re}_${as}_qde_char
         } 
      }
      
      foreach y in /*$basiccovs*/ $dumlist_meff ${grouplist_meff_`cb'} $conlist_meff age {
         if "`y'"~="agesq" {
            local me `me' me_`y'
         }
      }
      
      *Keeping only the marginal effects and percentage variation
      keep country `me' `perc' 

      *Calculating the variance of the generated average probabilities and marginal effects
      *One has to start from the 2nd line, since the first line gives the magnitudes
      *evaluated at the estimated coefficients
      local ncpp=${nsim}+1
      foreach vr in `me' `perc' {
         qui su `vr' in 2/`ncpp' 
         qui gen double se_`vr'=sqrt(r(Var))
      }

      *Calculating the point estimate as the mean of the effects across draws
      *We don't change the point estimate of the ownership proportions and their differences
      *since we use the actual ones
      foreach vr in `me' `perc' {
         if "`vr'"~="pr_${re}_${as}_`cb'_own" & "`vr'"~="dpus_${re}_${as}_`cb'_tot" & ///
            "`vr'"~="dpde_${re}_${as}_`cb'_tot" & "`vr'"~="pr_${re}_${as}_`cb'_cus" &  ///
            "`vr'"~="dpus_${re}_${as}_`cb'_coeff" & "`vr'"~="dpus_${re}_${as}_`cb'_char" & ///
            "`vr'"~="pr_${re}_${as}_`cb'_cde" & "`vr'"~="dpde_${re}_${as}_`cb'_coeff" & ///
            "`vr'"~="dpde_${re}_${as}_`cb'_char" {
            qui su `vr' in 2/`ncpp', meanonly 
            qui replace `vr' = r(mean) in 1
         } 
      }
  
      *Saving the dataset with all the probabilities and marginal effects
      qui save "${projf}/Results/SHARE/${re}/${as}/boot_alldraweff_${re}_${as}_`cb'.dta", replace

      *Saving the dataset with only the mean magnitudes and their standard errors
      qui keep in 1
      qui save "${projf}/Results/SHARE/${re}/${as}/boot_meaneff_${re}_${as}_`cb'.dta", replace
 
      *Saving matrices of marginal effects and percentage variations as a text file

      foreach rs in $prfx  {
         *Number of elements
         local n_`rs': word count ``rs''
         matrix olout=J(`n_`rs'',6,.)

         forvalues k=1/`n_`rs'' {
            local h: word `k' of ``rs''
            matrix olout[`k',1] = `h'[1]
            matrix olout[`k',2] = se_`h'[1]
            matrix olout[`k',3] = olout[`k',1]/olout[`k',2]
            matrix olout[`k',4] = 2*ttail(nobs_`cb' - ${nc},abs(olout[`k',3]))
            matrix olout[`k',5] = olout[`k',1] - invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
            matrix olout[`k',6] = olout[`k',1] + invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
         } 
    
         *Labels of matrix of magnitudes
         *Column labels
         local colmat `rs' se tstat pvalue ci_low ci_high
         matrix colnames olout = `colmat' 
         *Row labels
         matrix rownames olout = ``rs''
         *Displaying and saving the matrix of results
         mat list olout,format(%15.5f)
         mat2txt, mat(olout) ///
            saving("${projf}/Results/SHARE/${re}/${as}/boot_mefftab_${re}_${as}_`cb'.txt") replace format(%15.5f) ///
            title("Meff for ${as} in ${fnre}, using `cb' quartiles")

         erase "${projf}/Results/SHARE/${re}/${as}/temp2.dta"

      }  /* End of loop over magnitudes */        
   } /* End of loop over quartile specifications */
   
   erase "${projf}/Results/SHARE/${re}/${as}/temp.dta"         
   }  /*End of loop over countries */    
}  /*End of loop over assets */    


display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"
display "log file for ${as} in ${fnre} closed on $S_DATE at $S_TIME"
cap log close
   
   
  


